Simple gene signature to assess murine fibroblast polarization

We provide an original multi-stage approach identifying a gene signature to assess murine fibroblast polarization. Prototypic polarizations (inflammatory/fibrotic) were induced by seeded mouse embryonic fibroblasts (MEFs) with TNFα or TGFß1, respectively. The transcriptomic and proteomic profiles were obtained by RNA microarray and LC-MS/MS. Gene Ontology and pathways analysis were performed among the differentially expressed genes (DEGs) and proteins (DEPs). Balb/c mice underwent daily intradermal injections of HOCl (or PBS) as an experimental murine model of inflammation-mediated fibrosis in a time-dependent manner. As results, 1456 and 2215 DEGs, and 289 and 233 DEPs were respectively found in MEFs in response to TNFα or TGFß1, respectively. Among the most significant pathways, we combined 26 representative genes to encompass the proinflammatory and profibrotic polarizations of fibroblasts. Based on principal component analysis, this signature deciphered baseline state, proinflammatory polarization, and profibrotic polarization as accurately as RNA microarray and LC-MS/MS did. Then, we assessed the gene signature on dermal fibroblasts isolated from the experimental murine model. We observed a proinflammatory polarization at day 7, and a mixture of a proinflammatory and profibrotic polarizations at day 42 in line with histological findings. Our approach provides a small-size and convenient gene signature to assess murine fibroblast polarization.

www.nature.com/scientificreports/ chemo-attraction of immune cells and degradation of ECM. In response to tissue injury, some fibroblasts also acquire a myofibroblast phenotype, characterized by the expression of contractile proteins, and enhance ECM biosynthesis 2 . This matrix producer specialization is tightly controlled in normal conditions while in fibrotic diseases fibroblasts remain activated. The TGFβ/Smad pathway plays a central role in the profibrotic polarization of fibroblasts 4,5 , which is characterized by an excessive biosynthesis of ECM, an over-crosslinking, and a decrease of ECM degradation processes disrupting the tissue architecture and leading to fibrosis [6][7][8] . In an experimental murine model of inflammation-mediated fibrosis induced in a time-dependent manner, proinflammatory and profibrotic fibroblasts coexisted. The proinflammatory polarization predominated at early stage, while the profibrotic polarization predominated at the late stage 9,10 . The profibrotic polarization of fibroblasts is usually based on the evaluation of the gene expression of Col1a1 and Acta2, as markers of myofibroblasts, a central subpopulation of fibroblasts involved in ECM-overproduction 4,11 . Yet, none is specific on its own. Col1a1 is both upregulated in the inflammatory stage of the healing process and in fibrosis 9,10,12 . Several recent studies suggest that ECM overproduction is not limited to Acta2 + myofibroblasts and that Acta2 should not be used as a single marker of myofibroblastic transformation 9,10,13,14 .
Interpretation based on a limited number of genes can therefore be over-simplistic and lead to a misinterpretation of the polarization of fibroblasts. Deeper approaches like single-cell RNA sequencing have highlighted the heterogeneity of fibroblasts in models of wound healing and of fibrosis identifying new potential markers of fibroblast polarizations 12,14,15 . For example, it has been reported that other isoforms than Col1a1 were more expressed especially within ECM-producing fibroblasts (i.e. Col13a +) and inflammatory fibroblasts (i.e. Col5a3 + , Col14a +) 14 . Yet, these omics tools are time-and cost-consuming making them not suitable to explore the polarization of fibroblasts in routine.
Defining the polarization of fibroblasts and its evolution during time is very important to better understand the kinetic and pathophysiology of inflammation and fibrosis and to assess the impact of drugs, especially in animal models. To be more complete than overly simplistic single markers and more accessible than full (single cell) transcriptomics, we designed this study to find an accessible and accurate gene signature reflecting the polarization of fibroblasts. Firstly, we induced in vitro proinflammatory or profibrotic polarizations by treating Mouse Embryonic Fibroblasts (MEFs) with proinflammatory or profibrotic cytokines. Secondly, we selected polarization-representative genes using a combined analysis of the transcriptome assessed by RNA microarrays and of the proteome assessed by mass spectrometry (LC-MS/MS). Finally, we applied this gene signature at the inflammatory stage and at the fibrotic stage on dermal mouse primary fibroblasts (dMPF) isolated from an experimental murine model of inflammation-mediated fibrosis in a time-dependent manner [16][17][18][19] .

Results
Mouse embryonic fibroblasts treated by TNFα or TGFß1 are representative of proinflammatory or profibrotic polarizations of fibroblasts. First, we assessed the expression of Col1a1, Acta2, Fn1, Mmp3, Il6 in MEFs after a stimulation with proinflammatory (TNFα) or profibrotic (TGFß1) cytokines at 6-h, 12-h, 18-h, and 24-h. We observed that the gene expression profile associated with "profibrotic effects" was more marked after a 24-h stimulation than after a shorter stimulation (upregulation of Col1a1, down regulation of Mmp3 and Il6). Moreover, the distinction between TNFα-and TGFβ1-treated MEFs was more pronounced after a 24-h stimulation, especially with respect to Mmp3, Fn1 and Col1a1 (See Supplementary Fig. S1). For these reasons, we investigated the transcript expression of MEFs by RNA microarrays after a 24-h stimulation. With an absolute fold-change of 1.5, 1456 and 2215 distinct transcripts were significantly deregulated (q < 0.05) in MEFs treated by TNFα or TGFß1 compared to controls, respectively. Principal component analysis (PCA) showed that TNFα and TGFß1 induced two different gene expression profiles (See Supplementary Fig. S2). Based on the hierarchical clustering, we then focused analyses on clusters of interest of Differentially Expressed Genes (DEGs): (i) a cluster named TGFß1_up g of 412 DEGs preferentially overexpressed in MEFs treated by TGFß1, and (ii) a cluster named TNFα_up g of 242 DEGs preferentially overexpressed in MEFs treated by TNFα (See Fig. 1a). In these clusters of interest, a clustering enrichment analysis was performed to take account of redundancies between GO, KEGG and REACTOME terms. The analysis on TGFß1_up g showed an enrichment in terms associated to the ECM organization (i.e., GO:0048729: tissue morphogenesis; R-MMU-1474244: extracellular matrix organization), while the analysis on TNFα_up g showed an enrichment in terms associated to chemoattraction signaling (i.e., GO:0050900: leukocyte migration), immune response (i.e., GO:0050778: positive regulation of immune response), and inflammation response (i.e., GO:0006954: inflammation response) (See Fig. 1b- We next investigated the proteome using LC-MS/MS after a 72-h stimulation. Among the expressed proteins, 289 and 243 were significantly deregulated in MEFs (DEPs, Differentially Expressed Proteins) treated by TNFα or TGFß1 compared to controls, respectively (See Fig. 2a; Supplementary Fig. S4a-b). The Spearman correlation coefficients between the transcriptomic and proteomic expression were high (r = 0.79-0.83) (See Fig. 2b-c). PCA showed that these two cytokines induced two distinct protein expression profiles (See Supplementary Fig. S4c-e). As for RNA microarrays, we focused on: a cluster named TGFß1_up prot of 57 DEPs preferentially over-expressed in MEFs treated by TGFß1 and a cluster named TNFα_up prot of 27 DEPs preferentially over-expressed in MEFs treated by TNFα. MEFs treated by TGFß1 were enriched in terms associated with TGFß1 signaling (i.e., R-MMU-445144: signal transduction by L1) and ECM-modeling (i.e., mmu05205: proteoglycan in cancer), while MEFs treated by TNFα were enriched in pathways related to acute phase response (i.e., GO:0042743: hydrogen peroxide metabolic process) (See Fig. 2d www.nature.com/scientificreports/ or TGFß1 were representative of proinflammatory or profibrotic polarizations, we designed a gene signature to screen the activation of proinflammatory and profibrotic pathways. As LC-MS/MS did not allow the identification of smaller molecules such as cytokines, which are important mediators in the induction of these profiles, we focused on the process of selecting candidate genes based on transcriptome data, while including in our selection process the proteins identified by LC-MS/MS except for the smaller molecules. We selected several pathways involved in the acute phase response (GO:0006954) and in the organization of ECM (R-MMU-1474244) with a particularly focus on the crosslinking of collagen fibrils (mmu2243919), collagen formation (R-MMU-1474290), ECM-receptor interaction (R-MMU-04512), acute phase response signaling (IPA), inhibition of Matrix Metalloproteases (IPA), proto-myofibroblastic transformation (IPA), degradation of ECM (R-MMU-1474228.1), chemokine signaling (mmu04062) and leukocyte transendothelial migration (mmu04670) (See Fig. 3; Supplementary Fig. S5). We applied the following criteria to select candidate genes: significance of q-value (vs. untreated MEFs), pathways membership (representativeness of meaningful biological processes), importance of the fold change (vs. untreated MEFs), distinctiveness of the 2 prototypic polarizations (i.e., up-regulated in one and down-regulated in the other), and if applicable the identification of the corresponding protein in the LC-MS/MS approach. The conventional markers (Acta2, Fn1, Col1a1, Tgfß1, IL6, Timp1) were systematically included. The full gene set was composed of: C3, Ccl2, Cxcl1, Cxcl5, Icam1 (chemoattraction); Dcn, Il6, Tnfα, Sod2 (acute phase response signaling pathway); Mmp3 (ECM degradation); Tfpi2, Timp1 (regulation of ECM degradation); Col1a1, Col5a1, Col5a3, Col7a1 (collagens); Fn1, Jag1, Loxl3, Pcolce2, TGFß1, Dcn (ECM organization); Itga5, Itgb3 (ECM binding); Acta2, Itga11, Tagln (myofibroblastic transformation) (See Fig. 4a). The Spearman coefficients between the gene set expression assessed by RNA microarray and by RT-qPCR were high (range from 0.90 to 0.98) suggesting a very good correlation between the two methods. PCA showed that the gene signature assessed by RT-qPCR was able to separate the two polarizations of fibroblasts as accurately as RNA microarray did (See Fig. 4b-d; Supplementary Figs. S6-S7). We further analyzed the protein expression of some of these candidate markers either by ELISA assays in supernatant samples (MMP3, CXCL1) or by Western  www.nature.com/scientificreports/ 7 (mean ± SD = 37.1% ± 4.4 vs. 25.9% ± 5.9, p = 0.004), increasing progressively over time in HOCl mice (at day 42: 56% ± 4.9 vs. 27.7 ± 5.9 in PBS) (See Fig. 5c). At day 42 in HOCl mice, we observed a loss of normal skin architecture in accordance with skin fibrosis at the late stage (See Fig. 5d-e). The percentage of leukocytes among all skin cells was higher in HOCl mice compared to PBS mice on day 7 (mean ± SD = 33.7% ± 9.5 vs 13.0% ± 4.7, respectively, p = 0.002), then gradually decreased at day 42 (14.3% ± 7.1 vs 8.3% ± 4.0, respectively, p = 0.14) in line with an early inflammatory stage (n = 4-5/group). We investigated the gene signature expression of dMPFs isolated from skin biopsy samples. At day 7, we observed a proinflammatory polarization of dMPFs in d7 HOCL mice (n = 4, 1 experiment) compared to d7 PBS mice (n = 4, 1 experiment) characterized by: (i) the over-expression of the genes related to chemoattraction, acute phase response signaling and ECM degradation; and (ii) the down-expression of the genes related to ECM biosynthesis and myofibroblast transformation. Tnfα was under-expressed in d7 HOCL mice compared to d7 PBS mice. At day 42, we observed a mixture of a proinflammatory and profibrotic polarizations of MPFs in d42 HOCL mice compared to d42 PBS mice characterized by: (i) few inflammatory features, remaining lower than those observed in d7 HOCL mice, excepted a marked up-regulation of Tnfα expression (p < 0.001); (ii) the overexpression of genes related to ECM-biosynthesis and ECM-remodeling; and (iii) the down-expression of genes related to ECM-degradation (See Fig. 6, Supplementary Fig. S11). . Gene expression of pathways related to the inflammatory response and to ECM organization in TGFß1-and TNFα-treated MEFs. Gene expression of pathways in TGFß1-and TNFα-treated MEFs expressed as log2 (Fold change vs. controls) assessed by RNA microarrays. Gene lists come from molecular signatures database (MSigDB) (http:// softw are. broad insti tute. org/ gsea/ msigdb/ index. jsp) and Ingenuity Pathway Analysis (Qiagen, Inc) (https:// digit alins ights. qiagen. com). Only significant genes are presented in heatmaps (q-value < 0.05) generated using Prism (https:// www. graph pad. com/ scien tific-softw are/ prism/). Candidate genes and classical genes are respectively marked with green or black mark. www.nature.com/scientificreports/

Discussion
Nowadays, transcript gene expression of fibroblast polarization is mostly based either on the expression analysis of a few genes with the bias of missing a significant effect due to inaccurate pre-selected targets, or on more indepth analyzes such as RNA microarrays or RNA sequencing which take time and resources. This study provides an original process to propose a convenient and feasible gene signature to analyze fibroblast polarization. First, we selected candidate genes by following an original multi-stage pathological process. To induce the profibrotic fibroblast in vitro, we treated MEFs with TGFβ1 conventionally used as a positive control for fibrosis 4,20 . To induce the proinflammatory profile in vitro, we treated MEFs with TNFα which is produced by macrophages during the inflammatory phase of wound healing 21 . The transcript gene expression and the gene ontology analysis associated with the proinflammatory and profibrotic MEFs induced in vitro were consistent with previous data on profibrotic and proinflammatory polarizations of fibroblasts 4,10,22 . We selected the most representative pathways, chose representative genes of each pathway and combined them to illustrate the proinflammatory and profibrotic polarizations. Second, we assessed the gene signature on dMPF isolated from an experimental murine model of inflammation-mediated fibrosis in a time-dependent manner. Maria et al. showed that daily injections of HOCl triggered inflammation and led to cellular polymorphous infiltrates (inflammatory stage), following by a progressive development of established fibrosis, less inflammatory, and made of disorganized collagen fibers destructuring all skin layers 19 . At the inflammatory stage, we observed a proinflammatory polarization of the dMPF, while we observed a mixture of proinflammatory and profibrotic polarization at the fibrotic stage in consistence with the histological findings and the pathogenesis of fibrosis in this model [16][17][18][19] . www.nature.com/scientificreports/ The gene signature was composed both of conventionally-used genes (sometimes with a low power of discrimination such as Acta2) and innovative genes selected based on the importance of the fold change, the membership in meaningful biological processes and on the differential deregulation between the proinflammatory or profibrotic polarizations.
Several studies have shown that the ECM overproduction is not limited to Acta2 + fibroblasts and that Acta2 should not be used as a single marker of myofibroblastic transformation 9,10,13,14 . That is the reason why we selected Tagln 9,12,14 and Itga11 to increase the detection of myofibroblastic transformation [23][24][25] . In our study, we observed an increased expression of Acta2 at the level of gene and protein expressions in the two prototypical polarizations, which highlighted that an interpretation based on a few markers could lead to a misinterpretation www.nature.com/scientificreports/ of fibroblast polarization. It was therefore important to add other features to assess the profibrotic polarization, such as matrix crosslinking, matrix degradation, and matrix microenvironment. Matrix crosslinking regulate matrix stiffness [26][27][28] . Lysyl oxidases (LOXs) and lysyl oxidase-like 1-4 (LOXL 1-4) are a group of enzymes that catalyzes cross-linking of collagens, thereby rendering these matrix proteins unable to be degraded. They have been shown to be positively correlated with fibrosis in numerous organs including kidney, heart, liver, and lung 10,11,27,29,30 . Additionally, recent studies have shown that LOX(L) inhibition can lead to reduction in activated fibroblasts, in collagen content and in crosslink formation, making Loxl3 a marker of interest 31,32 . Degradation and inhibition of the degradation of the ECM should be also considered to explore fibroblast polarization. We selected Mmp3 and 2 inhibitors of MMPs to assess this balance. A potent inhibitor of metalloproteases is Tfpi2 which cleaves the serine proteases involved in the cleavage of biologically active pro-MMPs. Timp1 is one of the other major inhibitors of metalloproteases overexpressed in fibrotic diseases 33 . The ECM is a very complex three-dimensional structure that is not only limited to a mechanical support role but whose elements are also involved in various signaling pathways and cellular processes, including inflammation, proliferation, apoptosis, and angiogenesis. Among these elements, proteoglycans critically provide structure to the ECM, but also regulate many signaling pathways, such as growth factors, cytokines, chemokines, and matrix-cell contacts 6,26,34 .
Decorin is a crucial proteoglycan, also called "the guardian from the matrix", that "decorates" collagen fibrils and regulates the cell cycle by trapping TGFß 34,35 . Decorin has also the ability to act as a pan-receptor inhibitor of tyrosin kinases, including epidermal growth factor receptor, vascular endothelial growth factor receptor 2, mesenchymal-epithelial transition factor receptor; thus restraining angiogenesis 26,36 . The selection of decorin brings new insights in the assessment of fibroblast polarization. In a work on primary murine fibroblasts isolated from healthy lung or experimental pulmonary fibrosis, Akamatsu et al. showed that: (i) "inflammatory" fibroblasts significantly overexpressed genes associated with collagens and chemokines compared to quiescent fibroblasts; (ii) myofibroblasts overexpressed genes associated with ECM-synthesis and crosslinking, and underexpressed genes associated with chemokines compared to quiescent fibroblasts 10 . We selected several relevant inflammatory markers which have already been reported www.nature.com/scientificreports/ at the inflammatory phase of wound healing and of bleomycin-induced alveolitis to assess the inflammatory features associated with proinflammatory polarization 10,21 .
Previous studies on murine fibroblasts showed Col1a1 and Col5a1 were overexpressed both in inflammatory fibroblasts and in myofibroblasts 12,14 , which is consistent with our results. Some collagens seem to be preferentially expressed in myofibroblasts (Col7a1) or inflammatory fibroblasts (Col5a3, Col14a1) 14 . We have therefore selected several types of collagens to better describe this heterogeneity. In a single-cell analysis performed at day 21 in a mouse model of pulmonary fibrosis induced by bleomycin, Xie et al. identified two subclusters among Col1a1 + fibroblasts: Col13a1 matrix fibroblasts with profibrotic features and Col14a1 matrix fibroblasts with proinflammatory features such as up-regulation of Col5a3, Cxcl12, Mmp3 and Dcn 14 . These results are very consistent with features of proinflammatory and profibrotic polarizations delineated by our gene signature, while the gene expression of Col14a1 and Col13a1 were not significant in this study. Overall, our results suggest that this gene signature could be meaningful, and the small size of the gene list makes it a convenient way to precisely assess fibroblast polarizations, for example to screen drugs. An analysis based solely on the conventionally-used genes would not have illustrated neither the importance of the proinflammatory features of fibroblasts at d7 (Il6, but also chemoattraction factors), nor the profibrotic multidimensional features of fibroblasts at d42 (activating of Notch pathway, proto-myofibroblastic transformation, ECM biosynthesis, negative regulation of ECM degradation). Our work has some limitations. We developed the gene signature from gene expression of MEFs. Our application on dMPFs confirms its applicability, but further work will be necessary to assess its interest in other organs and in tumoral microenvironment. Indeed, we did not assess the gene signature to explore cancer associated fibroblast (CAF) polarization. However, single cell RNA sequencing of several cancer types demonstrated the presence of transcriptionally distinct CAF population 37,38 . Across various cancer types, myofibroblastic (αSMAhigh) CAFs are associated with a matrix-producing contractile phenotype, whereas inflammatory (αSMA-low) CAFs are generally specialized in inflammatory cytokine and chemokine secretion [37][38][39][40] . Given the similarities between, on the one hand, the profibrotic fibroblast with the myofibroblastic CAF, and one the other hand, the proinflammatory fibroblast with the inflammatory CAF, it could be relevant to explore the CAF polarization in the tumoral microenvironment using the gene signature.
We provide here an original detailed multi-stage process to develop a gene signature to analyze fibroblast polarization, including its confirmation in an experimental model of inflammatory fibrosis. This gene signature will allow to be more precise than the usual single markers currently used and more feasible and convenient that the complex full transcriptomic approach.

Mouse embryonic fibroblasts. Cells. MEFs were purchased from American Type Culture Collection
(Balb/c/3T3 clone A31, ATCC, MD, USA) and were used between 8 and 15th passages. Experimental murine model of inflammation-mediated fibrosis. Animals. Six-week-old female BALB/c mice (Janvier Labs, France) were used in all experiments. Animals were housed in a specific pathogenfree facility, within autoclaved ventilated cages with sterile food and water ad libitum, under constant room temperature and with 12-h day-night cycles. This study was carried out in accordance with the local and national guidelines (directive #68/609 CEE). The protocol was approved by the Regional Ethics Committee on Animal Experimentation.
Experimental procedure. Experimental model was induced by daily intradermal injections of 300 μl of an HOCl-generating solution into the shaved backs of mice, using a 27-gauge needle and a 1-ml syringe, as previously described. The HOCl-generating solution was extemporaneously prepared by adding NaClO solution (9.6% as active chlorine) to a 100 mM KH2PO4 solution (pH 6.2). The NaClO amount was determined by measuring the optical density (OD) of the solution at 280 nm, and then adjusted to obtain an OD between 0.7 and 0.9. Control mice received injections of 300 μl of sterilized phosphate-buffered saline (PBS).
Sample collection. Mice were sacrificed by cervical dislocation under deep CO2 anesthesia at either day 7 (inflammatory stage) or day 42 (fibrotic stage) after the first injection. Skin samples were collected at the time of euthanasia near the injection site.
Dermal mouse primary fibroblasts culture. To isolate dMPFs, the skin was harvested, minced into small pieces, and incubated in PBS containing dispase II (5 mg/ml) overnight at 4 °C. After removing epidermis, the dermis was incubated in 1X PBS containing collagenase IV (4 mg/ml) for 1 h. Cells were pelleted by centrifugation, washed twice, and grown in DMEM containing 10% FBS. Cells were used at the second passage.
Histological evaluation. Skin samples embedded in paraffin were sliced into serial 4-μm sections. Dermal thickness at the injection site was assessed by performing a hematoxylin eosin (HE) staining and measuring the distance between the epidermal-dermal junction and the dermal-subcutaneous fat junction at a 20-fold mag- www.nature.com/scientificreports/ nification using the ImageJ software (U.S. National Institute of Health) (https:// imagej. nih. gov/ ij/ index. html). Twenty random measurements per section were performed by two blinded investigators and averaged for each section.
Measurement of collagen deposition. Collagen deposition in the skin was evaluated using Picrosirius red staining. The color was prepared by a 5-min bath of hematoxylin diluted 1/3 in water (Bio Optical). A 60 min bath of Picrosirius Red, which specifically dyes collagen fibers, was then carried out, obtained by mixing Direct red 80 (Sigma Aldrich) diluted to 1% in saturated picric acid (Sigma Aldrich). The sections were mounted between slide and coverslips and acquired under an automated microscope as described above, then analyzed using the ImageJ software (U.S. National Institute of Health) (https:// imagej. nih. gov/ ij/ index. html). A semi-quantitative analysis of the collagen deposition was performed using a method of Color deconvolution and the area occupied by collagen was quantified (expressed in %).
Leukocytes count in skin sample. After collection, the skin biopsies were digested with type II dispase (1 mg/ mL to 1.88 units/mL) and type IV collagenase (0.5 mg/mL to 215 units/mL) for 2h30 at 37 °C. The suspension was then centrifuged at 300 g for 10 min. RNA microarray data analysis. We analyzed RNA microarray data using Linear Models for Microarray Analysis (LIMMA) package in R software (version R 3.6.0) (https:// www.r-proje ct. org) to identify genes that were expressed differentially between datasets. The threshold for statistical significance was set to q-value (Benjamini-Hochberg false discovery rate (FDR) at < 0.05. Normalization with Z-score was determined, and hierarchical clustering was constructed with Euclidean distance. Ontology annotation and canonical pathway enrichment. To analyze Gene Oncology (GO), REACTOM, and KEGG pathways, we subjected differentially expressed genes (DEGs) (or proteins, DEPs) with cut-off q-value at < 0.05 to the online tool Metascape (http:// metas cape. org/ gp/ index. html#/ main/ step1) 43 . All significant enriched terms (p-value < 0.05) were hierarchically clustering into a tree based on Kappa-statistical similarities (threshold > 0.3) among their gene membership. The clusters identified have been named according to the membership. Among the terms within each cluster, some of them (q-value < 0.05) were selected for charting. Molecular signatures database (MSigDB) (http:// softw are. broad insti tute. org/ gsea/ msigdb/ index. jsp) and Ingenuity Pathway Analysis (Qiagen, Inc) (https:// digit alins ights. qiagen. com) were used to identify candidate genes in selected gene set.

Statistical analysis.
Values were reported as median ± interquartile range or mean ± SD according to the normality of distribution checked graphically and using the Shapiro-Wilk test. Comparisons of groups were performed using the analysis of variance or the Kruskal-Wallis test. In the case of significance results, the pairwise comparisons were performed. Comparisons were performed using Student t test for gaussian continuous variables and using Mann Whitney U test for non-gaussian continuous variables. Correlation was performed using Spearman correlation. All statistical tests were performed at 2-tailed α level of 0.05. All statistical analyses were performed on Prism (version 9.0.2 GraphPad Software) (https:// www. graph pad. com/ scien tific-softw are/ prism/). Venn's diagrams were created using Venny 2.1 (BioinfoGP-Csic) (https:// bioin fogp. cnb. csic. es/ tools/ venny).
Ethics approval and consent to participate. This study was carried out in accordance with the local and national guidelines (directive #68/609 CEE). The protocol was approved by the Regional Ethics Committee on Animal Experimentation CEEA 75 (Comité d' éthique sur l' expérimentation animale) under the reference number: APAFIS#19603-2020061914271271 v6.
Additional statement. The study was carried out in compliance with the ARRIVE guidelines.

Data availability
Transcriptomics datasets related to this article can be found at https:// www. www.nature.com/scientificreports/